home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / arith / idiv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  10.1 KB  |  512 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  idiv.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16.  
  17. /*    Umstellung auf PLACEadd etc., RD, 24.9.91    */
  18. /*    Normalisieren fuer die Division durch Shifts.    */
  19. /*    31.8.92, RD, Fehler in Idiv korrigiert: bb, one nicht
  20.         zerstoert (Hinweis von H.B. Eggenstein vom 10.8.92)
  21. */
  22.  
  23. #include <LEDA/impl/iint.h>
  24. #include <LEDA/impl/iloc.h>
  25.  
  26. static void normdiv(u, v, q, ul, vl)
  27.     register pPLACE u, v;
  28.     pPLACE q;
  29.     int ul;
  30.     register int vl;
  31. /* u, v sind verschiedene Integer.vec der Laengen ul, vl >=2, die fuer 
  32.    die Division u/v normiert sind. Das heisst auch, dass 
  33.    v[vl-1] > u[ul-1] ist. u[ul-1]=0 ist zugelassen. v ist nicht 0. 
  34.    Zurueckgegeben wird der Rest in u und der Quotient in q, die Laenge
  35.    von q muss groesser oder gleich ul-vl+1 sein. */
  36. {   register int j;
  37.     for (j=ul-1; j>=vl; j--) {
  38.     PLACE qdach, h, l, carry, help;
  39.     if (u[j]==v[vl-1]) {
  40.         qdach=RADIXMINUSONE;
  41.         /* (h, l) = (v[vl-2] - u[j-1])*RADIX - u[j-2] - 1 */
  42.         carry=PLACEsub(0, u[j-2], 1, &l);
  43.         carry=PLACEsub(v[vl-2], u[j-1], carry, &h);
  44.         if (!carry) {
  45.         /* (h, l) -= (v[vl-1]*RADIX + v[vl-2]) */
  46.         carry=PLACEsub(l, v[vl-2], carry, &l);
  47.         carry=PLACEsub(h, v[vl-1], carry, &h);
  48.         while (!carry) {
  49.             qdach-=1;
  50.             carry=PLACEsub(l, v[vl-2], carry, &l);
  51.             carry=PLACEsub(h, v[vl-1], carry, &h);
  52.         }
  53.         }
  54.     } else {
  55.         help=PLACEdiv(u[j], u[j-1], v[vl-1], &qdach);
  56.         /* (h, l) = v[vl-2]*qdach - help*RADIX - u[j-2] - 1; */
  57.         h=PLACEmul(qdach, v[vl-2], 0, &l);
  58.         carry=PLACEsub(l, u[j-2], 1, &l);
  59.         carry=PLACEsub(h, help, carry, &h);
  60.         while (!carry) {
  61.         qdach-=1;
  62.             carry=PLACEsub(l, v[vl-2], carry, &l);
  63.             carry=PLACEsub(h, v[vl-1], carry, &h);
  64.         }
  65.     }
  66.     /* Nun ist qdach bestimmt, und hoechstens 1 zu gross */
  67.     carry=vecmulsub(u+j-vl, v, qdach, vl);
  68.     if (carry) {        /* Addiere v zurueck */
  69.         qdach-=1;
  70.         carry=vecaddto(u+j-vl, v, vl);
  71.         PLACEadd(u[j], 0, carry, &u[j]);    /* Ignoriere carry */
  72.     }
  73.     q[j-vl]=qdach;
  74.     }        /* j-Schleife */
  75.     return;
  76. }        /* normdiv */
  77.  
  78. static PLACE norm_vecsli(u, a, toshift, count)
  79.     PLACE *u, *a;
  80.     int toshift, count;
  81. /* 0<=toshift<=LOGPLACE-1; */
  82. {    PLACE accu, help;
  83.     int i, bleft = LOGPLACE - toshift;
  84.  
  85.     if ( ! toshift ) {
  86.         for (i=0; i<count; i++)
  87.             u[i]=a[i];
  88.         return 0;
  89.     }
  90.     accu = 0;
  91.     for ( i=0; i<count; i++ ) {
  92.         help = a[i];
  93.         accu |= (help<<toshift);
  94.         u[i]=accu;
  95.         accu=help>>bleft;
  96.     }
  97.     return accu;
  98. }        /* norm_vecsli */
  99.  
  100. static void norm_vecsri(a, u, toshift, count)
  101.     PLACE *a, *u;
  102.     int toshift, count;
  103. /* 0<=toshift<=LOGPLACE-1; */
  104. {    PLACE accu, help;
  105.     int i, bleft=LOGPLACE-toshift;
  106.  
  107.     if ( ! toshift ) {
  108.         for (i=0; i<count; i++)
  109.             a[i]=u[i];
  110.         return;
  111.     }
  112.     accu=u[0];
  113.     accu>>=toshift;
  114.     for ( i=1; i<count; i++ ) {
  115.         help=u[i];
  116.         accu|=(help<<bleft);
  117.         a[i-1]=accu;
  118.         accu=help>>toshift;
  119.     }
  120.     a[count-1]=accu;
  121. }        /* norm_vecsri  */
  122.  
  123.  
  124. /*******************************************************/
  125.  
  126. #ifdef __GNUC__
  127. PLACE uIdiasP(Integer *a, PLACE b)
  128. #else
  129. PLACE uIdiasP(a, b)
  130.     register pInteger a;
  131.     register PLACE b;
  132. #endif
  133. {    register PLACE rem;
  134.     if (!a->length)
  135.         return 0;
  136.     rem=vecdiv(a->vec, a->vec, b, a->length);
  137.     if (!a->vec[a->length-1])
  138.         a->length--;
  139.     return rem;
  140. }    /* uIdiasP */
  141.  
  142. #ifdef __GNUC__
  143. PLACE IdiasP(Integer *a, PLACE b)
  144. #else
  145. PLACE IdiasP(a, b)
  146.     register pInteger a;
  147.     register PLACE b;
  148. #endif
  149. {    register PLACE rem;
  150.     Integer one;
  151.     if (!a->length)
  152.         return 0;
  153.     rem=vecdiv(a->vec, a->vec, b, a->length);
  154.     if (!a->vec[a->length-1])
  155.         a->length--;
  156.     if (a->sign==PLUS)
  157.         return rem;
  158.     if (!rem)
  159.         return rem;
  160.     cIasint(&one, 1);
  161.     ImiasI(a, &one);
  162.     dI(&one);
  163.     return b-rem;
  164. }    /* IdiasP */
  165.  
  166. #ifdef __GNUC__
  167. PLACE uIasIdiP(Integer *q, const Integer *a, PLACE b)
  168. #else
  169. PLACE uIasIdiP(q, a, b)
  170.     register pInteger q;
  171.     register const Integer *a;
  172.     register PLACE b;
  173. #endif
  174. {    register PLACE rem;
  175.     register int nl;
  176.     if (q==a)
  177.         return uIdiasP(q, b);
  178.     if (!a->length) {
  179.         Iasint(q, 0);
  180.         return 0;
  181.     }
  182.     nl=a->length;
  183.     if (nl>q->maxlength) {
  184.         delvec(q->vec, q->maxlength);
  185.         q->maxlength=nl;
  186.         q->vec=newvec(&q->maxlength);
  187.     }
  188.     rem=vecdiv(q->vec, a->vec, b, a->length);
  189.     if (q->vec[a->length-1])
  190.         q->length=a->length;
  191.     else
  192.         q->length=a->length-1;
  193.     q->sign=PLUS;
  194.     return rem;
  195. }    /* uIasIdiP */
  196.  
  197. #ifdef __GNUC__
  198. PLACE IasIdiP(Integer *q, const Integer *a, PLACE b)
  199. #else
  200. PLACE IasIdiP(q, a, b)
  201.     register pInteger q;
  202.     register const Integer *a;
  203.     register PLACE b;
  204. #endif
  205. {    register PLACE rem;
  206.     register int nl;
  207.     Integer one;
  208.     if (q==a)
  209.         return IdiasP(q, b);
  210.     if (!a->length) {
  211.         Iasint(q, 0);
  212.         return 0;
  213.     }
  214.     nl=a->length;
  215.     if (nl>q->maxlength) {
  216.         delvec(q->vec, q->maxlength);
  217.         q->maxlength=nl;
  218.         q->vec=newvec(&q->maxlength);
  219.     }
  220.     rem=vecdiv(q->vec, a->vec, b, a->length);
  221.     if (q->vec[a->length-1])
  222.         q->length=a->length;
  223.     else
  224.         q->length=a->length-1;
  225.     if (a->sign==PLUS) {
  226.         q->sign=PLUS;
  227.         return rem;
  228.     }
  229.     q->sign=MINUS;
  230.     if (!rem)
  231.         return rem;
  232.     cIasint(&one, 1);
  233.     ImiasI(q, &one);
  234.     dI(&one);
  235.     return b-rem;
  236. }    /* IasIdiP */
  237.  
  238. /*************************************************/
  239.  
  240. void uIdiv(q, r, a, b)
  241.     pInteger q, r;
  242.     const Integer *a, *b;
  243. /* Division mit Rest, a=bq+r. */
  244. {   register pPLACE u, v;
  245.     register int ul, vl;
  246.     int unl, vnl, toshift;
  247.     register int m;
  248.     register int i;
  249.     PLACE help, carry;
  250.  
  251.     vl=b->length;
  252.     m=a->length-vl;
  253.     if (m<0) {
  254.     IasI(r, a);
  255.     r->sign=PLUS;
  256.     Iasint(q, 0);
  257.     return;
  258.     }
  259.     if (vl<=1) {
  260.     register PLACE rem;
  261.     if (!vl)
  262.         Ierror("uIdiv: division by zero");
  263.     rem=uIasIdiP(q, a, b->vec[0]);
  264.     *(r->vec)=rem;
  265.     if (rem)
  266.         r->length=1;
  267.     else
  268.         r->length=0;
  269.     r->sign=PLUS;
  270.     return;
  271.     }
  272.     /* Hilfsvariablen bereitstellen */
  273.     ul=a->length+1;
  274.     unl=ul;
  275.     u=newvec(&unl);
  276.     vnl=vl;
  277.     v=newvec(&vnl);
  278.     /* a, b normalisieren */
  279.     help=b->vec[vl-1];
  280.     toshift=0;
  281.     while ( !( help >> (LOGPLACE-1)) ) {
  282.     help<<=1;
  283.     toshift++;
  284.     }
  285.     u[ul-1]=norm_vecsli(u, a->vec, toshift, ul-1);
  286.     norm_vecsli(v, b->vec, toshift, vl);
  287.     /* eigentliche Division */
  288.     if (m+1>q->maxlength) {
  289.     delvec(q->vec, q->maxlength);
  290.     q->maxlength=m+1;
  291.     q->vec=newvec(&q->maxlength);
  292.     }
  293.     normdiv(u, v, q->vec, ul, vl);
  294.     /* Rest zurueckgewinnen */
  295.     if (vl>r->maxlength) {
  296.     delvec(r->vec, r->maxlength);
  297.     r->maxlength=vl;
  298.     r->vec=newvec(&r->maxlength);
  299.     }
  300.     norm_vecsri(r->vec, u, toshift, vl);
  301.     delvec(u, unl);
  302.     i=vl-1; 
  303.     u=r->vec;
  304.     while(!u[i] && i>=0)
  305.     i--;
  306.     r->length=i+1;
  307.     r->sign=PLUS;
  308.     /* Quotient auf Standardform */
  309.     if (q->vec[m])
  310.     q->length=m+1;
  311.     else
  312.     q->length=m;
  313.     q->sign=PLUS;
  314.     delvec(v, vnl);
  315. }        /* uIdiv */
  316.  
  317. /*************************************************/
  318.  
  319. void Idiv(q, r, a, b)
  320.     pInteger q, r;
  321.     const Integer *a, *b;
  322. /* Division mit Rest, a=bq+r. */
  323. {   register pPLACE u, v;
  324.     register int ul, vl;
  325.     int unl, vnl, toshift;
  326.     register int m;
  327.     register int i;
  328.     PLACE help, carry;
  329.     BOOLEAN usebb=FALSE;
  330.     Integer bb, one;
  331.  
  332.     vl=b->length;
  333.     m=a->length-vl;
  334.     if (m<0) { /* dann sind a, b verschiedene Variablen */
  335.     if (a->sign==PLUS) {
  336.         IasI(r, a);
  337.         Iasint(q, 0);
  338.         return;
  339.     }
  340.     if (b->sign==PLUS) {
  341.         IasIplI(r, a, b);
  342.         Iasint(q, -1);
  343.         return;
  344.     } else {
  345.         IasImiI(r, a, b);
  346.         Iasint(q, 1);
  347.         return;
  348.     }    
  349.     }        /* m<0 */
  350.     /* Sonderfall: Divisor einstellig */
  351.     if (vl<=1) {
  352.     register PLACE rem;
  353.     if (!vl)
  354.         Ierror("Idiv: Division by zero");
  355.     rem=IasIdiP(q, a, b->vec[0]);
  356.     *(r->vec)=rem;
  357.     if (rem)
  358.         r->length=1;
  359.     else
  360.         r->length=0;
  361.     r->sign=PLUS;
  362.     q->sign^=b->sign;
  363.     return;
  364.     }
  365.     /* Hilfsvariablen bereitstellen */
  366.     ul=a->length+1;
  367.     unl=ul;
  368.     u=newvec(&unl);
  369.     vnl=vl;
  370.     v=newvec(&vnl);
  371.     /* a zu u, b zu v normalisieren */
  372.     help=b->vec[vl-1];
  373.     toshift=0;
  374.     while ( !( help >> (LOGPLACE-1)) ) {
  375.     help<<=1;
  376.     toshift++;
  377.     }
  378.     u[ul-1]=norm_vecsli(u, a->vec, toshift, ul-1);
  379.     norm_vecsli(v, b->vec, toshift, vl);
  380.     /* eigentliche Division */
  381.     if (a->sign==MINUS) {
  382.     if ((b==r)||(b==q)) {
  383.         usebb=TRUE;
  384.         cIasI(&bb, b);
  385.     }
  386.     }
  387.     if (m+1>q->maxlength) {
  388.     delvec(q->vec, q->maxlength);
  389.     q->maxlength=m+1;
  390.     q->vec=newvec(&q->maxlength);
  391.     }
  392.     normdiv(u, v, q->vec, ul, vl);
  393.     /* Rest zurueckgewinnen */
  394.     if (vl>r->maxlength) {
  395.     delvec(r->vec, r->maxlength);
  396.     r->maxlength=vl;
  397.     r->vec=newvec(&r->maxlength);
  398.     }
  399.     norm_vecsri(r->vec, u, toshift, vl);
  400.     delvec(u, unl);
  401.     i=vl-1; 
  402.     u=r->vec;
  403.     while(!u[i] && i>=0)
  404.     i--;
  405.     r->length=i+1;
  406.     /* Quotient auf Standardform */
  407.     if (q->vec[m])
  408.     q->length=m+1;
  409.     else
  410.     q->length=m;
  411.     delvec(v, vnl);
  412.     /* Rest positiv, auf a==bq+r normalisieren. */
  413.     if (a->sign==PLUS) {
  414.     if (!(q->length))
  415.         q->sign=PLUS;
  416.     else
  417.         q->sign=b->sign;
  418.         r->sign=PLUS;
  419.         if (usebb)
  420.          {dI(&bb);}
  421.     return;
  422.     }
  423.     if (!r->length) {
  424.     if (!q->length)
  425.         q->sign=PLUS;
  426.     else
  427.         q->sign=a->sign^b->sign;
  428.     r->sign=PLUS;
  429.         if (usebb)
  430.          {dI(&bb);}
  431.     return;
  432.     }
  433.     cIasint(&one, 1);
  434.     if (!usebb) {
  435.         if (b->sign==PLUS) {
  436.         r->sign=MINUS;
  437.         IplasI(r, b);
  438.         q->sign=PLUS;
  439.         IplasI(q, &one);
  440.         q->sign=MINUS;
  441.             dI(&one);
  442.         return;
  443.     } else {
  444.         r->sign=MINUS;
  445.         ImiasI(r, b);
  446.         q->sign=PLUS;
  447.         IplasI(q, &one);
  448.             dI(&one);
  449.         return;
  450.     }
  451.     } else {
  452.         if (bb.sign==PLUS) {
  453.         r->sign=MINUS;
  454.         IplasI(r, &bb);
  455.         q->sign=PLUS;
  456.         IplasI(q, &one);
  457.         q->sign=MINUS;
  458.             dI(&bb); 
  459.             dI(&one);
  460.         return;
  461.     } else {
  462.         r->sign=MINUS;
  463.         ImiasI(r, &bb);
  464.         q->sign=PLUS;
  465.         IplasI(q, &one);
  466.             dI(&bb); 
  467.             dI(&one);
  468.         return;
  469. }   }   }    /* Idiv */
  470.  
  471. /*****************************************/
  472.  
  473. void IasIdiI(q, a, b)
  474.     pInteger q;
  475.     const Integer *a, *b;
  476. /* Division, Quotient */
  477. {    Integer r;
  478.     cI(&r);
  479.     Idiv(q, &r, a, b);
  480.     dI(&r);
  481. }
  482.  
  483. void IdiasI(q, b)
  484.     pInteger q;
  485.     const Integer *b;
  486. /* Division, Quotient */
  487. {    Integer r;
  488.     cI(&r);
  489.     Idiv(q, &r, q, b);
  490.     dI(&r);
  491. }
  492.  
  493. void IasIreI(r, a, b)
  494.     pInteger r;
  495.     const Integer *a, *b;
  496. /* Division, Rest */
  497. {    Integer q;
  498.     cI(&q);
  499.     Idiv(&q, r, a, b);
  500.     dI(&q);
  501. }
  502.  
  503. void IreasI(r, b)
  504.     pInteger r;
  505.     const Integer *b;
  506. /* Division, Rest */
  507. {    Integer q;
  508.     cI(&q);
  509.     Idiv(&q, r, r, b);
  510.     dI(&q);
  511. }
  512.